Multiscaling at Point J: Jamming is a Critical Phenomenon 
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We analyze the jamming transition that occurs as a function of increasing packing density in 
a disordered two-dimensional assembly of disks at zero temperature for "Point J" of the recently 
proposed jamming phase diagram. We measure the total number of moving disks and the transverse 
length of the moving region, and find a power law divergence as the packing density increases toward 
a critical jamming density. This provides evidence that the T = jamming transition as a function 
of packing density is a second order phase transition. Additionally we find evidence for multiscaling, 
indicating the importance of long tails in the velocity fluctuations. 
PACS numbers: 45.70.-n, 64.60.Ht 



There has been a surge of activity in jamming phenom- 
ena for T = systems such as granular materials, foams, 
and colloids, where jamming is defined as the onset of a 
nonvanishing yield stress in a disordered state [1,2]. Liu 
and Nagel proposed a jamming phase diagram containing 
three axes: temperature T, the inverse packing fraction 
1/0, and shear a [1]. The system is jammed within a 
three-dimensional dome; above the jamming transition, 
the system behaves as a rigid solid. Recently O'Hern et 
al. [3] studied the area of the jamming phase diagram 
at T = a = as a function of packing density (f>, and 
showed that a well defined sharp jamming transition ap- 
pears at "Point J." Elsewhere on the jamming phase di- 
agram, the boundary between jammed and unjammed 
states is not sharp since its definition is sensitive to the 
experimental time scale. It was noted in Ref. [3] that, 
although Point J does not exist for liquids, the behavior 
near Point J could be relevant to the physics near the 
glass transition. A key question is whether Point J is a 
true continuous phase transition with a diverging spatial 
correlation length. Since the physics at Point J is non- 
thermal, any continuous phase transition behavior could 
be dominated by rare fluctuations. If this is the case, 
simple scaling close to the critical density may not occur 
even if there is an underlying phase transition. 

Here, we consider the motion of a single disk driven 
through an arrangement of disks as the jamming tran- 
sition is approached, a system which, as suggested in 
Ref. [3], should provide a direct test of the nature of 
the jamming transition. The reader is invited to try the 
following experiment: place a large collection of coins flat 
on a desk, so that they are almost touching. Then, push 
one coin and observe what other coins move. At low 
packing fractions 4>, the driven disk or coin occasionally 
contacts other disks, which can in turn contact still other 
disks, but the total number of disks moving is small. At a 
certain density 4> = (f> c , however, the entire system must 
move simultaneously as a unit, and jams. Force chains 
can be observed emanating from the driven disk [4] . 

In these jamming problems, the tail of the particle ve- 
locity distribution P(v) at small velocities controls the 



behavior of the system, since once the system enters a 
jammed state, with everything stopped in the thermody- 
namic limit, it can never exit the jammed state. The best 
method for characterizing the tail of P(v) is via multi- 
scaling, which we employ in this Letter. For fixed driving 
force, a small v implies that a large number of other par- 
ticles are dragged with the driven particle, but we can 
also characterize the number of particles moved by the 
driven particle by examining the force transmission. We 
will see that this measurement provides direct evidence 
for a diverging length scale. Our data also suggests that 
the jamming transition for very large systems coincides 
with the random close packed density. 

We simulate a binary mixture of two-dimensional disks 
with stiff spring repulsive interactions of radius ta and 
r B at T = (see Fig. 1). For all the densities we consider 
here, we find <C 1% overlap in the radii as the disks in- 
teract, indicating that the spring constant is sufficiently 
large to provide a good approximation to hard core disks. 
The bimodal disk distribution, with a size ratio of 1.4 : 1, 
is chosen to create a disordered arrangement and avoid 
formation of a regular lattice. We also performed sim- 
ulations with a size ratio of 1.7 : 1 with substantially 
identical results. We employ overdamped dynamics such 
that the velocity of each disk is proportional to the force 
acting on it. The equation of motion for the disks is 



r)Vi = Vfe(|ry| -r e}f )^- +f d . 

■ / ■ \ r ij 



(1) 



Here k 

[5], r i:j - 



200 is the strength of the stiff restoring spring 
Yi — Yj , and for all but one of the disks, the ex- 
ternal driving force = 0. The single driven disk (large 
dot in Fig. 1) has = 0.1. Disks only interact if they 
are separated by a distance smaller than the sum of their 
radii, |ry| < r e // = r^' +r^>, where r^' equals either ta 
or r B , depending on the given particle. We tested several 
different values of Lj and chose the value small enough to 
minimize overlap as mentioned above. We use periodic 
boundary conditions with a range of linear system sizes 
L = 24 to 60. For L — 60, the system contains N = 2600 
background disks at a density of <fi = 0.8395. 
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FIG. 1. Sample geometry for i = 60 and (a) <f> = 0.656, 
(b) 4> = 0.811, (c) (f> = 0.837. Large black dot: the probe disk 
(size exaggerated for clarity); gray dots: moving disks; small 
dots: non-moving disks. 

We prepare the system by randomly dropping nonover- 
lapping disks up to <j> w 0.6. To reach higher densi- 
ties, we add disks at randomly selected interstitial lo- 
cations, reduce the radii of all disks, and then increase 
the radii back to the initial values while allowing the sys- 
tem to evolve under a small temperature. This produces 
nonoverlapping configurations at densities up to 4> c . To 
reach the maximum possible density, v / 37r/6 w 0.907, re- 
quires phase separating the system. At the significantly 
lower (f> c fa 0.839 we find no phase separation. Since the 
system is not jammed below <p c , we believe that we suc- 
cessfully uniformly sampled all available states, and that 
phase separation does not occur until some 4> sep > <p c . 

To study the velocity curve, we drive a single large disk 
along a 45° angle over a distance of V2L/5. The time re- 
quired for this motion is much longer than the time scale 
of any brief initial transients in the velocity. We measure 
the effective velocity of the driven disk over the length 
of the simulation. For computational reasons, it is not 
possible to simulate a disk moving at an arbitrarily slow 
velocity through the system. In these cases, we declared 
that the system had jammed; this definition agreed well 
with the critical <f> c w 0.839 given by curve fitting below. 

In Fig. 1 we show snapshots of simulation results for 
different densities, with disks counted as moving and 
shaded in gray if they were connected via a force contact 
to the driven disk. At the lowest <j> = 0.656, in Fig. 1(a), 
the driven disk moves easily and interacts with only a 
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FIG. 2. Time series of velocity for L = 60 and <f> — (a) 
0.656, (b) 0.747, and (c) 0.837. 

few neighbors. In Fig. 1(b) for <f> = 0.811, close to jam- 
ming, a larger portion of the disks are moving, while in 
Fig. 1(c) at <j) = 0.837, the system is sufficiently close to 
4> c that the entire (finite-size) sample is moving. 

Figure 2 shows example time series for the velocity v 
parallel to the applied drive at <f> = 0.656, 0.747, and 
0.837. Not only does the average velocity decrease as <ft 
increases, but the time series becomes more intermittent. 
At <j> = 0.837, the velocity is usually very small, but there 
are occasional bursts of much higher velocity. 

We emphasize that for < <p c the threshold force re- 
quired to move the driven particle vanishes. Instead, v 
is proportional to the drive fd at any moment of time. 
This follows from dimensional grounds, since for a hard- 
core system, there is no physical parameter with units of 
force. This is very different from the behavior shown in 
a recent numerical study of single probe particles driven 
through a system with softer, long-range forces [7], where 
there is a non- vanishing, but finite, depinning force at all 
densities due to the long-range nature of the force. 

Distribution of Velocities — Newton's third law, com- 
bined with Eq. (1), leads to the result that r]J2i v i = 
Thus, if the driven particle is not in contact with any 
other particles, it moves at v = fd/V- If there are a total 
of n particles moving together, including the driven par- 
ticle, then they each move at v = fd/n. For <j> > <p c , all 
N particles in the system move together, and v is van- 
ishingly small in the thermodynamic limit. For 4> < <p c , 
the number of particles moving is finite in the thermody- 
namic limit and diverges as <f> approaches cf> c . 

One measurement which may show scaling as <j> ap- 
proaches <j) c is the average velocity v of the driven disk. 
However, v is a random, time-dependent variable, and for 
any density <fi < <f> c there is a nonvanishing probability of 
observing any given, arbitrarily small v. Since the other 
disks are distributed randomly, there is still some proba- 
bility that in any local region the density will exceed 4> c 
[8] (this argument is inspired by Lifshitz tails). If this 
region encloses n disks, the probability of finding such a 
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FIG. 3. Plot of M(q,4>) 1/q versus cj> for q = -1, 1, 2, and 
4 (from bottom to top). The curves with larger q are nois- 
ier. Inset: M{q,4>) x ^ q versus <j> c — 4> on a log-log scale for 
ij = —3, —1, 1, 3. A curve fit shows the power law scaling. 

region is exponentially small in n, but it is nonvanishing. 
When the driven disk hits such a region, it behaves as 
if it jams, and v slows to a value of order 1/n until the 
disk either escapes the region or pushes other disks into 
neighboring regions of lower density. Thus, with expo- 
nentially small probability in n, v ~ 0(l/n). Once the 
particle slows, however, it takes a long time for it to leave 
the region. Thus, one expects to see an intermittent v. 

To measure the intermittency, we use multifractal [9] 
scaling. For a given 0, let p(v)dv be the probability of 
measuring a given v at a single instant in time. Dchne 
the q-ih inverse moment by M(q) = J dvp(v)v~ q as the 
time average of the q-th power of the inverse velocity. We 
then define a set of multifractal exponents r(q) by 

M(q) = J dvp(v)v- q cx (<f> c - 4>)- T(q) . (2) 

If v were constant throughout the simulation for a given 
4>, we would have M(q) 1 / q — M(l). If instead v had rela- 
tive fluctuations of order unity about some characteristic 
velocity, vo((f>), then M(q) 1 ^ q would not equal M(l) in 
general, but r(q)/q would still be independent of q. We 
will instead find that r(q)/q depends on q. This implies 
that v does not simply fluctuate about a single vo, but is 
instead much more intermittent, with the particle some- 
times getting stuck for a long time at a slow velocity, 
then traveling much more rapidly. 

We compute the moments M and extract the expo- 
nents r(q) by fitting the moments to the form c * (4> c — 
0)- T (<?) The value of <j) c determined from these fits is 
independent of q to high accuracy for q > —2, and thus 
we can fix <p c ss 0.839 as the onset of jamming. In Fig. 3 
we plot M(q, (j)) 1 / q versus </> for L = 60 and various q, av- 
eraged over three realizations of the system. In the inset 
we show the log-log plot of the curves in the main panel 




FIG. 4. (a): Plot of r(q)/q against q. The q dependence 
of r(q)/q shows the existence of multiscaling in this system, 
(b): Plot of M(l) (crosses) and Af(3) 1/3 (circles) against 
M(-l)- 1 . 

with solid lines indicating power law fits. The curves 
not only fail to overlap, but also have different slopes, 
indicating the presence of multiscaling. In Fig. 4(a), we 
plot r(q)/q versus q. For large positive and negative q, 
the plot asymptotes. The asymptote at large positive q 
reflects the scaling of the typical (disregarding exponen- 
tially rare possibilities discussed above) slowest velocity 
of the system which goes to zero as (</> — </> c ) llm «^=° T ( q )/ q . 
Similarly, the asymptote at negative q reflects the typi- 
cal largest velocity. Since lim^-oo r(q)/q is very close to 
zero, the typical largest velocity is largely independent of 
(j>, as seen in Fig. 2. In fact, it is likely that at sufficiently 
negative q the exponent r(q)/q becomes zero. 

The error in r(q) due to statistical fluctuations is neg- 
ligible. If we compare the r(q) obtained by curve fitting 
M(q, 4>) averaged over all realizations to r(q) obtained by 
using only a single realization, the difference in exponents 
is of order 0.001 to 0.002. Instead, the major source of 
error is finite size effects. The straight lines in the inset 
to Fig. 3 are based on fitting over a certain range of <j). 
For <j) closer to (p c than the endpoint of the fitting lines, 
the statistical noise in the data increases significantly and 
some of the realizations jam while others do not. Thus, 
for finite N the jamming threshold is not well defined, 
and these <fi are so close to <p c that finite size effects may 
be important. The error bars in Fig. 4 are based on this 
consideration of finite size effects. The low end of the er- 
ror bars corresponds to the r(q) obtained by fitting only 
over the range shown in the inset to Fig. 3, while the 
high end corresponds to including all <j> < <p c . The low 
end tends to underestimate the difference in r(q)/q for 
different q. r(q)/q is clearly g-dependent since the differ- 
ence in r(q)/q between q = —1 and q = 2, for example, 
is definitely outside the error bars. In Fig. 4(b) we plot 
M{qY/ q as a function of M(-l)" 1 for q = 1,3 over the 
same range as in Fig. 3 to illustrate that the data obey ex- 
tended self-similarity [10]. The slopes on the log-log plot 
are 1.25 ± 0.15 and 1.45 ± 0.2 respectively; as expected, 
this is within error bars of the ratio of the exponents 
-t(1)/t(-1) and -(r(3)/3)/r(-l). It is difficult to es- 
tablish that these slopes are different from each other, 
but the slopes are definitely different from unity, which 
already indicates multiscaling. The fact that r(— 1) is 
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FIG. 5. Plot of n mov ing against (f>. Solid line: L = 60. 
Dashed line: L — 48. Dotted line: a power law fit to the 
L = 60 curve. Left inset: The L — 60 curve and the fit against 
4> c — 4>, showing scaling. Right inset: finite size scaling. 

less than one makes the size of the scaling regime seem 
small, since a wide scaling range in (j> leads to only a small 
range in M(— 1); however, the curve fit is very good over 
the available range. 

Figure 1 shows that as the jamming transition is ap- 
proached, the number of moving disks increases and there 
is a diverging length scale as the jamming is approached. 
In order to quantify this, we show in Fig. 5 the number 
of moving disks n mov ing vs <f> for systems of linear size 
L = 48 and L = 60. Here a clear divergence appears as 
the critical density is approached. The divergence is cut 
off when n mov in g equals the total number of disks. We 
have also considered smaller systems for different param- 
eters of drive and disk radii and again observe a diver- 
gence; however, these smaller systems give a much lower 
resolution and hence a larger error on the estimated ex- 
ponent. In Fig. 5 we fit a power law to the largest system 
with (<f) c — cf)) T , where r is between 1.2 and 1.46. We also 
measured the sum over moving disks of the squared dis- 
tance between the driven disk and the moving disk, a 
quantity n mov ing^ 2 , and find that this number diverges as 
(4> c — 0)~ CT , with a between 2.35 and 2.6. The number of 
moving disks cannot directly be compared to the velocity 
v of the driven disk, as the driven disk can push other 
disks normal to the drive so that n mov ing may be much 
larger than 1/v. We obtained the exponents r and a by 
a time average of the number of moving disks, and did 
not perform a multiscaling calculation, though we can 
extract a length from a simple scaling analysis. If the 
moving disks form a cluster with dimension d and length 
scale £ diverging as £ oc (<fr c — <j>)~", then r = vd and 
a = v(d + 2). Given the values of r, a, and our own 
qualitative observations, we find d = 2, so that v is be- 
tween 0.6 and 0.7. This provides a direct measurement 
of the diverging length scale, in reasonable agreement 
with the finite-size scaling result [3] v = 0.71 ± 0.12. We 



note that Fig. 5 in our case is completely consistent with 
finite-size scaling as shown in the right inset. Scale the 
y-axis by L? (the scaling for the total number of parti- 
cles in the system), and scale <f> c — <j> by L~ 2 l r . Then, 
the curves for different L automatically collapse within 
the scaling region, and since the divergence is cutoff at 
"-moving <x. L 2 for all L, these curves also collapse close to 
<f> c . Finally, our (f> c rts 0.839 is very close to the critical 
packing density found in [3]. 

Conclusion — We have studied a system of a T = 0, 
zero shear 2D disordered assembly of disks at densities 
below and up to Point J in the recently proposed jam- 
ming phase diagram. A single probe disk is pushed with 
a constant drive through the other disks. Upon increas- 
ing the packing density, we find a jamming transition 
associated with a power law divergence in the number of 
moving disks and a diverging spatial correlation length, 
indicating that Point J is a true continuous phase tran- 
sition. In addition, we show that the tails of the disk 
velocity distribution play an important role in this tran- 
sition, since due to the non-thermal nature of the system, 
once the particle stops moving, it cannot restart. Using 
the multifractal moments r(q), we chart the tails and 
show that the r(q) do not depend linearly on q; thus, 
the system exhibits multiscaling. We also find evidence 
for a diverging correlation length of the force contacts as 
the jamming transition is approached from below. Our 
results show that there is an underlying second order non- 
thermal phase transition at the jamming transition. We 
suggest that experimental work should consider both the 
velocity time series and the number of moving particles in 
the surrounding media, and should test for multiscaling 
behavior as the jamming transition is approached. Par- 
ticular experiments would include driving a single disk 
on a flat surface through a disordered assembly of other 
disks for increasing density, or driving individual colloids 
through glassy assemblies of other colloids. 
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